library(mosaic)
# install.packages("remotes") # avkommentera första gången
library(remotes)
#install_github("StatisticsSU/sda1paket") # avkommentera första gången
library(sda123)Statistik och dataanalys I, 15 hp
Datorlaboration 7 - Inferens för regression
Installera paket
Den här labben förutsätter att följande paket finns installerade:
mosaicsda123
sda123 är kursens egna R-paket och måste installeras från GitHub. Det görs genom att först installera paketet remotes och därefter installera kurspaketet:
Introduktion
I den här datorlabben kommer ni analysera olika datamaterial med regression. Det första datamaterialet om reklam är verkligt. I Uppgift 3 ska ni analysera datamaterial som är simulerade av mig, från olika populationsmodeller som jag har valt. Det blir ett intressant laboratorium, där ni kan utforska olika aspekter av regression i en kontrollerad miljö. Och även få veta den underliggande populationsmodellen efter att ni gjort labben! The truth will reveal itself!
Skapa mapp för labben
Skapa en mapp Lab7 i din kursmapp SDA1. Ladda ner Quarto-filen för denna lab genom att högerklicka här och välj ‘Spara länk’ eller något likande från menyn som dyker upp. Spara filen i din Lab7 mapp. Öppna Quarto-filen i RStudio och fortsätt med denna laboration direkt i Quarto-dokumentet, där du också fyller i svaren på dina laborationsövningar.
1. Enkel regression för reklamdata
Datamaterialet reklam.csv innehåller data på n=200 produkters försäljning (sales) i tusentals US dollar och hur mycket (i s k reklamenheter) produkten har marknadsförts i olika reklamkanaler: tv, radio och newspaper. Vi läser in data med kommandot:
reklam = read.csv(file = 'https://github.com/StatisticsSU/SDA1/raw/main/datorlab/lab7/reklam.csv')
head(reklam) tv radio newspaper sales
1 230.1 37.8 69.2 22.1
2 44.5 39.3 45.1 10.4
3 17.2 45.9 69.3 12.0
4 151.5 41.3 58.5 16.5
5 180.8 10.8 58.4 17.9
6 8.7 48.9 75.0 7.2Vi börjar med en enkel regressionsmodell med enbart tv som förklarande variabel. Vi utforskar data genom histogram för responsvariabeln sales och scatterplot mot tv.
par(mfrow = c(1,2))
hist(reklam$sales, col = "orange", main = "")
plot(sales ~ tv, data = reklam, pch = 19, cex = 0.7, col = "steelblue")
💪 Uppgift 1.1
Skatta den enkla linjära regressionsmodellen i R:
\[ \texttt{sales} = \beta_0 + \beta_1 \cdot \texttt{tv} + \varepsilon, \hspace{1cm} \varepsilon\overset{iid}{\sim}N(0,\sigma_{\varepsilon}) \]
och tolka skattningen av intercept och lutningskoefficient. Vilken skattning av \(\sigma_{\varepsilon}\) ger R, och hur tolkar du denna?
Den skattningen av intercept är ca 6.97 och lutningskoefficient är ca 0.06 och skattning av \(\sigma_{\varepsilon}\) är ca 2.30.
model1 = lm(sales ~ tv, data = reklam)
summary(model1)
Call:
lm(formula = sales ~ tv, data = reklam)
Residuals:
Min 1Q Median 3Q Max
-6.4438 -1.4857 0.0218 1.5042 5.6932
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.974821 0.322553 21.62 <2e-16 ***
tv 0.055465 0.001896 29.26 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.296 on 198 degrees of freedom
Multiple R-squared: 0.8122, Adjusted R-squared: 0.8112
F-statistic: 856.2 on 1 and 198 DF, p-value: < 2.2e-16💪 Uppgift 1.2
Testa om variabeln tv är en signifikant förklarande variabel på 5% signifikansnivå. Ställ upp noll- och alternativhypotes, teststatistiska med fördelning under nollhypotesen. Utför testet med hjälp av formelsamlingen och formulera din slutsats. Du får använda information från R, men du ska ställa upp formeln för teststatistikan och beräkna med den formeln.
model1 = lm(sales ~ tv, data = reklam)
confint(model1) 2.5 % 97.5 %
(Intercept) 6.33874038 7.61090260
tv 0.05172671 0.05920283summary(model1, conf_intervals = TRUE, anova = FALSE)
Call:
lm(formula = sales ~ tv, data = reklam)
Residuals:
Min 1Q Median 3Q Max
-6.4438 -1.4857 0.0218 1.5042 5.6932
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.974821 0.322553 21.62 <2e-16 ***
tv 0.055465 0.001896 29.26 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.296 on 198 degrees of freedom
Multiple R-squared: 0.8122, Adjusted R-squared: 0.8112
F-statistic: 856.2 on 1 and 198 DF, p-value: < 2.2e-16t_teststatistikan <- (0.06-0)/0.002
t_teststatistikan[1] 30n = nrow(reklam)
qt(0.975,n-2)[1] 1.972017t_teststatistikan > qt(0.975,n-2) # if TRUE -> reject H0: Beta1 = 0[1] TRUE💪 Uppgift 1.3
Visa hur R har beräknat p-värdet i regression-utskriften (från summary) genom att göra beräkningen själv med relevant s k p,d,q, eller r-funktion.
model1 = lm(sales ~ tv, data = reklam)
summary(model1)
Call:
lm(formula = sales ~ tv, data = reklam)
Residuals:
Min 1Q Median 3Q Max
-6.4438 -1.4857 0.0218 1.5042 5.6932
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.974821 0.322553 21.62 <2e-16 ***
tv 0.055465 0.001896 29.26 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.296 on 198 degrees of freedom
Multiple R-squared: 0.8122, Adjusted R-squared: 0.8112
F-statistic: 856.2 on 1 and 198 DF, p-value: < 2.2e-16p_value <- 2*pt(t_teststatistikan, n-2, lower.tail = F) #P(Tobs>=1.97)
p_value[1] 1.399538e-75pt(1.97, df=198)#sannolikhet[1] 0.9748837💪 Uppgift 1.4
Gör en residualanalys med funktionen reg_residuals i sda123-paketet, och kommentera om var och ett av de fyra grundläggande antaganden för populationsmodellen verkar uppfyllda.
Figur ett är qqplot, data ligger på theiretical linje och tolkar som normalfördelning.
library(sda123)
model1 = lm(sales ~ tv, data = reklam)
summary(model1)
Call:
lm(formula = sales ~ tv, data = reklam)
Residuals:
Min 1Q Median 3Q Max
-6.4438 -1.4857 0.0218 1.5042 5.6932
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.974821 0.322553 21.62 <2e-16 ***
tv 0.055465 0.001896 29.26 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.296 on 198 degrees of freedom
Multiple R-squared: 0.8122, Adjusted R-squared: 0.8112
F-statistic: 856.2 on 1 and 198 DF, p-value: < 2.2e-16reg_residuals(model1)
💪 Uppgift 1.5
Gör en prediktion med 95%-igt prediktionsintervall för sales vid tv=100. [tips: argumentet newdata i predict-funktionen måste vara en dataframe, inte en vektor. Se min kod för lifespan data.]
model1 = lm(sales ~ tv, data = reklam)
reg_summary(model1)
Analysis of variance - ANOVA
------------------------------------------------
df SS MS F Pr(>F)
Regr 1 4512.4 4512.4352 856.18 7.9279e-74
Error 198 1043.5 5.2704
Total 199 5556.0
Measures of model fit
------------------------------------------------
Root MSE R2 R2-adj
2.29575 0.81218 0.81123
Parameter estimates
------------------------------------------------
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.974821 0.3225535 21.624 5.0277e-54
tv 0.055465 0.0018956 29.260 7.9279e-74# Prediktionsintervall - sda123
reg_predict(sales ~ tv, data = reklam)`geom_smooth()` using formula = 'y ~ x'
qt(p=0.975, df=198) #sannolikhet[1] 1.972017# Anpassat värde/prediktion för Sverige:
predict(model1, newdata = data.frame(tv = 100),interval="prediction") fit lwr upr
1 12.5213 7.979338 17.063262. Multipel regression för reklamdata
Vi går nu vidare med en multipel regressionsmodell med alla tre förklarande variabler: tv, radio och newspaper. Först en titt på data:
par(mfrow = c(2,2))
hist(reklam$sales, col = "orange", main = "")
plot(sales ~ tv, data = reklam, pch = 19, cex = 0.7, col = "steelblue")
plot(sales ~ radio, data = reklam, pch = 19, cex = 0.7, col = "steelblue")
plot(sales ~ newspaper, data = reklam, pch = 19, cex = 0.7, col = "steelblue")
💪 Uppgift 2.1
Anpassa följande multipla regressionsmodell:
\[ \texttt{sales} = \beta_0 + \beta_1 \cdot \texttt{tv} + \beta_2 \cdot \texttt{radio} + \beta_3 \cdot \texttt{newspaper} + \varepsilon, \hspace{1cm} \varepsilon\overset{iid}{\sim}N(0,\sigma_{\varepsilon}) \]
och tolka skattningen av regressionskoefficienten för tv.
model2 <- lm(sales ~ tv + radio + newspaper, data = reklam)
summary(model2)
Call:
lm(formula = sales ~ tv + radio + newspaper, data = reklam)
Residuals:
Min 1Q Median 3Q Max
-7.3034 -0.8244 -0.0008 0.8976 3.7473
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.6251241 0.3075012 15.041 <2e-16 ***
tv 0.0544458 0.0013752 39.592 <2e-16 ***
radio 0.1070012 0.0084896 12.604 <2e-16 ***
newspaper 0.0003357 0.0057881 0.058 0.954
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.662 on 196 degrees of freedom
Multiple R-squared: 0.9026, Adjusted R-squared: 0.9011
F-statistic: 605.4 on 3 and 196 DF, p-value: < 2.2e-16💪 Uppgift 2.2
Testa om variabeln tv är en signifikant förklarande variabel på 5% signifikansnivå. Ställ upp noll- och alternativhypotes, teststatistiska med fördelning under nollhypotesen. Utför testet och formulera din slutsats. Du får använda information från R, men du ska ställa upp formeln för teststatistikan och beräkna den.
Kommentera kortfattat utifrån utskriften om radio och newspaper är signifikanta på 5% signifikansnivå.
Radio är signifikanta på 5% signifikansnivå utan newspaper inte är signifikanta på 5% signifikansnivå.
model2 <- lm(sales ~ tv + radio + newspaper, data = reklam)
confint(model2) 2.5 % 97.5 %
(Intercept) 4.01868836 5.23155980
tv 0.05173372 0.05715784
radio 0.09025861 0.12374384
newspaper -0.01107921 0.01175052summary(model2)
Call:
lm(formula = sales ~ tv + radio + newspaper, data = reklam)
Residuals:
Min 1Q Median 3Q Max
-7.3034 -0.8244 -0.0008 0.8976 3.7473
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.6251241 0.3075012 15.041 <2e-16 ***
tv 0.0544458 0.0013752 39.592 <2e-16 ***
radio 0.1070012 0.0084896 12.604 <2e-16 ***
newspaper 0.0003357 0.0057881 0.058 0.954
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.662 on 196 degrees of freedom
Multiple R-squared: 0.9026, Adjusted R-squared: 0.9011
F-statistic: 605.4 on 3 and 196 DF, p-value: < 2.2e-16# H0: beta_tv = 0, Ha: beta_tv != 0
t_teststatistikan <- (0.0544458-0)/0.0013752
t_teststatistikan[1] 39.59119n = nrow(reklam)
qt(0.975,n-2)[1] 1.972017t_teststatistikan > qt(0.975,n-2) # if TRUE -> reject H0: Beta1 = 0[1] TRUE💪 Uppgift 2.3
Vilken modell, den enkla regressionen eller den multipla, föredrar du? Motivera med relevant R² värde och genom att göra 5-fold korsvalidering med funktionen reg_crossval i sda123-paketet.
RMSE_model2 = reg_crossval(sales ~ tv + radio + newspaper, data = reklam, nfolds = 5)
RMSE_model2 [1] 1.6773743. Residualanalys på simulerade data
Jag har simulerat fem olika datamaterial simdata1-5 från en linjär regressionsmodell:
\[ \texttt{y} = \beta_0 + \beta_1 \cdot \texttt{X1} + \beta_2 \cdot \texttt{X2} + \varepsilon \]
där i vissa datamaterial har feltermer \(\varepsilon\) som inte uppfyller alla av de fyra antagandena, och eventuellt kan det också finnas ett icke-linjärt samband mellan någon förklarande variabel och y. Er uppgift är att utföra residualanalys med funktionen reg_residuals i sda123-paketet för att försöka upptäcka vilka antaganden som inte är uppfyllda i respektive datamaterial (Uppgift 3.1-3.5 nedan). Ni bör också titta på plot(simdata1) för varje datamaterial för att upptäcka icke-linjära samband.
Här är de fem datamaterialen:
simdata1 = read.csv(file = 'https://github.com/StatisticsSU/SDA1/raw/main/datorlab/lab7/simdata1.csv')
simdata2 = read.csv(file = 'https://github.com/StatisticsSU/SDA1/raw/main/datorlab/lab7/simdata2.csv')
simdata3 = read.csv(file = 'https://github.com/StatisticsSU/SDA1/raw/main/datorlab/lab7/simdata3.csv')
simdata4 = read.csv(file = 'https://github.com/StatisticsSU/SDA1/raw/main/datorlab/lab7/simdata4.csv')
simdata5 = read.csv(file = 'https://github.com/StatisticsSU/SDA1/raw/main/datorlab/lab7/simdata5.csv')💪 Uppgift 3.1 - simdata1
head(simdata1) X y X1 X2
1 1 1.2856459 0.30665435 0.66858877
2 2 0.4728535 -0.59227811 0.29043835
3 3 0.4576525 -1.48905857 0.21114563
4 4 -0.1214482 -0.64102846 -0.55785480
5 5 -3.2053896 -1.20232219 -0.09721979
6 6 -0.3920066 0.04123726 1.71397704m1 = lm(y ~ X1 + X2, data = simdata1)
reg_residuals(m1)
summary(m1)
Call:
lm(formula = y ~ X1 + X2, data = simdata1)
Residuals:
Min 1Q Median 3Q Max
-12.3951 -0.5728 0.1052 0.8056 5.6914
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.9851 0.1101 8.945 2.74e-16 ***
X1 1.0369 0.1035 10.014 < 2e-16 ***
X2 -0.1385 0.1067 -1.298 0.196
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.541 on 197 degrees of freedom
Multiple R-squared: 0.3513, Adjusted R-squared: 0.3447
F-statistic: 53.35 on 2 and 197 DF, p-value: < 2.2e-16plot(simdata1)
💪 Uppgift 3.2 - simdata2
m2 = lm(y ~ X1 + X2, data = simdata2)
reg_residuals(m2)
summary(m2)
Call:
lm(formula = y ~ X1 + X2, data = simdata2)
Residuals:
Min 1Q Median 3Q Max
-3.5471 -0.5977 0.0276 0.6443 4.0607
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.03782 0.07845 13.230 <2e-16 ***
X1 0.99503 0.08783 11.329 <2e-16 ***
X2 0.04638 0.07719 0.601 0.549
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.108 on 197 degrees of freedom
Multiple R-squared: 0.3952, Adjusted R-squared: 0.3891
F-statistic: 64.38 on 2 and 197 DF, p-value: < 2.2e-16plot(simdata2)
💪 Uppgift 3.3 - simdata3
m3 = lm(y ~ X1 + X2, data = simdata3)
reg_residuals(m3)
summary(m3)
Call:
lm(formula = y ~ X1 + X2, data = simdata3)
Residuals:
Min 1Q Median 3Q Max
-2.67225 -0.53425 0.09314 0.55420 2.02347
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.03070 0.06272 16.434 <2e-16 ***
X1 0.92495 0.06195 14.931 <2e-16 ***
X2 0.09723 0.06693 1.453 0.148
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.8854 on 197 degrees of freedom
Multiple R-squared: 0.5352, Adjusted R-squared: 0.5305
F-statistic: 113.4 on 2 and 197 DF, p-value: < 2.2e-16plot(simdata3)
💪 Uppgift 3.4 - simdata4
m4 = lm(y ~ X1 + X2, data = simdata4)
reg_residuals(m4)
summary(m4)
Call:
lm(formula = y ~ X1 + X2, data = simdata4)
Residuals:
Min 1Q Median 3Q Max
-5.3459 -1.4186 -0.0414 1.3921 5.0268
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.68543 0.14816 4.626 6.73e-06 ***
X1 0.83993 0.13810 6.082 6.07e-09 ***
X2 0.02523 0.13708 0.184 0.854
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.077 on 197 degrees of freedom
Multiple R-squared: 0.1582, Adjusted R-squared: 0.1496
F-statistic: 18.51 on 2 and 197 DF, p-value: 4.302e-08plot(simdata4)
💪 Uppgift 3.5 - simdata5
m5 = lm(y ~ X1 + X2, data = simdata5)
reg_residuals(m5)
summary(m5)
Call:
lm(formula = y ~ X1 + X2, data = simdata5)
Residuals:
Min 1Q Median 3Q Max
-0.8418 -0.2863 0.0024 0.2579 1.0555
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.35010 0.02720 12.873 <2e-16 ***
X1 0.02583 0.04673 0.553 0.581
X2 -0.03729 0.02755 -1.354 0.177
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.3824 on 197 degrees of freedom
Multiple R-squared: 0.01134, Adjusted R-squared: 0.001303
F-statistic: 1.13 on 2 and 197 DF, p-value: 0.3252plot(simdata5)